R Notebook

# Import packages
library(tidyverse)
library(lubridate)
library(RQuantLib)
library(plotly)

Step 0: Import data and reference indices; remove superfluous variables to create datasets for agency correction

# Import needed indices
setwd('Datasets/Indices')
park_ids <- read_csv('Park and Agency Codes Index.csv') # For Step 2
entrance_classifications <- read_csv('2018 Entrance Usage Classifications Index.csv') # For Step 3 & 5
bus_vehicle_multipliers <- read_csv('Bus and Vehicle Multipliers Index.csv') # For Step 4

Run the following import and export once, prior to requesting Qualtrics accuracy audits at the end of the season

# Import 2018 Qualtrics datasets - see Qualtrics Info document for login information
# Note that for import by code below to work, all datasets must be csvs, saved with original Qualtrics names in 'Datasets' folder (i.e. download as csv from Qualtrics and do not change the name)
setwd('.')
setwd('Datasets/Qualtrics')
parks_18 <- list.files(pattern = '*.csv')

make_names <- function(name_vector) {
 
  name_vector_df <- as.data.frame(name_vector)
 
 new_names <- name_vector_df %>%
   rename(Agency_name = 1) %>%
   mutate(Agency_name = str_replace_all(Agency_name, c('\\+-+.+?v' = '', '\\+' = '_')))

 new_names_list <- as.list(new_names)
 
 return(new_names_list$Agency_name)
}

list2env(map(set_names(parks_18, make_names(parks_18)), read_csv), envir = .GlobalEnv)
## <environment: R_GlobalEnv>

all_agencies_list <- list(Anoka_County = Anoka_County, Bloomington = Bloomington, Carver_County = Carver_County, Dakota_County = Dakota_County, Minneapolis_Park_and_Recreation_Board = Minneapolis_Park_and_Recreation_Board, Ramsey_County = Ramsey_County, Saint_Paul = Saint_Paul, Scott_County = Scott_County, Three_Rivers_Park_District = Three_Rivers_Park_District, Washington_County = Washington_County)

# Remove unnecessary variables to clean datasets prior to sending to agencies
# 
# variable_paring <- function(df) {
#   df %>%
#     slice(-2) %>%
#     select(18:33)
# }
# 
# parks_pared <- map(all_agencies_list, variable_paring)

# Write out datasets without extra variables to folder on N drive, from which they can be sent for correction

# parks_pared %>%
#   names(.) %>%
#   map(~ write_csv(parks_pared[[.]], paste0("Datasets/To Send for Accuracy Audit/", ., ".csv")))

Step 1: Tidy agencies’ final data (with corrections), bind together

#setwd('.')
#setwd('Datasets/Audited')
#parks_18 <- list.files(pattern = '*.csv')

all_agencies_list <- list(Anoka_County, Bloomington, Carver_County, Dakota_County, Minneapolis_Park_and_Recreation_Board, Ramsey_County, Saint_Paul, Scott_County, Three_Rivers_Park_District, Washington_County)

park_tidy <- function(df) {
  df %>%
    slice(-1:-2) %>%
  rename(Counter_initials = 18,
         Agency_name = 19,
         Park_trail = 20,
         Entrance = 21,
         Month = 22,
         Day = 23,
         Year = 24,
         Timeslot = 25,
         Scheduled_count = 26,
         Scheduled_date = 27,
         Bikers = 28,
         Peds_skaters = 29,
         Vehicles = 30,
         Buses = 31,
         Horse_riders = 32,
         Boaters = 33) %>%
  select(8, 9, 18:33) %>%
  mutate(Start_time = ifelse(Timeslot != '10:00-12:00 PM', str_replace(Timeslot, '\\-.+? ', ''), '10:00AM')) %>%
  separate(Start_time, into = c('Start', 'AM_PM'), sep = -2) %>%
  unite(Start_time, Start, AM_PM, sep = ' ') %>%
  mutate(Zero = 0,
         Start_time_2 = Start_time) %>%
  unite(Start_alt, Zero, Start_time_2, sep = '') %>%
  mutate(Start_time_fmt = ifelse(Start_time != '10:00 AM' & Start_time != '12:00 PM', Start_alt, Start_time)) %>%
  select(-Start_alt)
  
}

agencies_tidied <- map(all_agencies_list, park_tidy)

agencies_bound <- do.call(bind_rows, agencies_tidied)

Step 2: Create a day type variable (weekend/holiday or weekday); create a variable with the total number of days for that day type in the season

Make sure, in this step, to change the year in the first two lines to the pertinent year - this has to be done manually because while usually use estimates are done the year after summer counts are complete, it could be that they are done in the same year as the counts in years to come

# Calculate days in the summer season for the year
memorial_labor <- as.data.frame(RQuantLib::getHolidayList(calendar = 'UnitedStates', from = as_date(ymd('20180501')), to = as_date(ymd('20180930')), includeWeekends = FALSE))

summer_days <- memorial_labor %>%
  rename(Date = 1) %>%
  slice(-2) %>% # remove Fourth of July
  mutate(row_id = row_number()) %>%
  mutate(Holiday = ifelse(row_id == 1, 'Memorial_Day', 'Labor_Day')) %>%
  select(-row_id) %>%
  spread(Holiday, value = Date) %>%
  mutate(Memorial_Day_Saturday = as_date(ymd(Memorial_Day))- days(2)) %>% # Summer season begins the Saturday before Memorial Day
  mutate(Season_period = as.period(interval(start = as_date(ymd(Memorial_Day_Saturday)), end = as_date(ymd(Labor_Day))), unit = 'day')) %>%
  separate(Season_period, into = c('Season_days', "HMS_in_season"), sep = 'd') %>%
  mutate(Half_of_summer = (as.numeric(Season_days)/2)) %>%
  mutate(Midpoint_of_summer = as_date(ymd(Memorial_Day_Saturday)) + days(Half_of_summer))

# Calculate number of weekend/holiday days and number of weekdays in summer season - this 
day_type_counts <- as.data.frame(seq(as.Date(summer_days$Memorial_Day_Saturday), as.Date(summer_days$Labor_Day), "days")) %>%
  rename(Date = 1) %>%
  mutate(Holiday = RQuantLib::isHoliday(calendar = "UnitedStates", dates = Date)) %>%
  mutate(Day_type = ifelse(Holiday == 'TRUE', 'Weekends_holidays', 'Weekdays')) %>%
  group_by(Day_type) %>%
  count() %>%
  spread(Day_type, value = n)

# Assign which half of the summer the day falls in (this is used, in part, ); assign day-type variable; assign number of days of that day-type in the season (this number is used as a multiplier to scale counts up)
step_II <- agencies_bound %>%
  mutate(Year = ifelse(is.na(Year), '2018', Year)) %>% # Some agencies don't enter the year; fill in for them
  mutate(Month_unite = Month,
         Day_unite = Day,
         Year_unite = Year) %>%
  unite(Date, Month_unite, Day_unite, Year_unite, sep = " ") %>%
  mutate(Date = as_date(mdy(Date))) %>%
  mutate(Half_of_summer = ifelse(Date <= summer_days$Midpoint_of_summer, 1, 2)) %>% 
  unite(Date_year, Month, Day, Year, sep = '-') %>%
  mutate(Date_object = as_date(mdy(Date_year))) %>%
  mutate(Holiday = RQuantLib::isHoliday(calendar = "UnitedStates", dates = Date_object)) %>%
  mutate(Day_type = ifelse(Holiday == 'TRUE', 'Weekend_holiday', 'Weekday')) %>%
  mutate(Total_days_by_type = ifelse(Day_type == 'Weekend_holiday', day_type_counts$Weekends_holidays, day_type_counts$Weekdays)) %>%
  select(-Holiday) %>%
  mutate(Park_trail = ifelse(Park_trail == 'Minnesota River Greenway', 'Minnesota River Greenway Regional Trail', Park_trail)) %>%
  mutate(Agency_name = ifelse(Park_trail == 'Battle Creek Regional Park', 'Ramsey County',
                              ifelse(Park_trail == 'Mississippi Gorge Regional Park', 'Minneapolis Park and Recreation Board',
                                     ifelse(Park_trail == 'Hidden Falls-Crosby Farm Regional Park', 'Saint Paul',
                                            ifelse(Park_trail == 'Scott County Regional Trail', 'Scott County', Agency_name))))) %>% # These agency selection for these trails from the Qualtrics form are incorrect
  mutate(Start_time = as.character(Start_time_fmt)) %>%
  select(-Start_time_fmt)

Step 2a: Add numeric Park Identifier to dataset (used in step 3 to add park entrance usage classifications)

step_IIa <- left_join(step_II, park_ids, by = c('Park_trail' = 'Park_name', 'Agency_name' = 'Agency_name'))

# Ensure all data joined properly

step_IIa_unjoined <- step_IIa %>%
  filter(is.na(Park_id) | is.na(Agency_id))

Step 2a check: Check the number of observations in the step_IIa dataset. If the number of observations increased, some of the counts have been duplicated upon join. Use the following code to discern why and how to ameliorate the faux pas.

# Which agency or agencies have duplicated counts?

# step_II_agency_counts <- step_II %>%
#   group_by(Agency_name) %>%
#   count()
# 
# step_IIa_agency_counts <- step_IIa %>%
#   group_by(Agency_name) %>%
#   count()
# 
# step_IIa_compare_agency_counts <- full_join(step_II_agency_counts, step_IIa_agency_counts, by = c('Agency_name'))
# 
# # In this case, only Anoka.  Which Anoka parks or trails were duplicated?
# 
# step_II_park_counts <- step_II %>%
#   filter(Agency_name == 'Anoka County') %>%
#   group_by(Park_trail) %>%
#   count()
# 
# step_IIa_park_counts <- step_IIa %>%
#   filter(Agency_name == 'Anoka County') %>%
#   group_by(Park_trail) %>%
#   count()
# 
# step_IIa_compare_park_counts <- full_join(step_II_park_counts, step_IIa_park_counts, by = c('Park_trail'))
# 
# step_IIa_compare_counts %>%
#   filter(n.x != n.y)
# # In this case, Rum River Central Regional Park is responsible for the duplication.  Rum River CRP was given two park attributions for Anoka County, by mistake.  Remove the '109' attribution from the index.
# 
# step_IIa_unjoined <- step_IIa %>%
#   filter(is.na(Park_id) | is.na(Agency_id)) %>%
#   select(Park_trail, Agency_name) %>%
#   unique()

Step 3: Add Park Entrance Usage Classifications (High, Medium, Low)

# Note that Trout Brook's RT Entrance 4 was fixed manually; if the 2018 Entrance Usage Classifications Index is re-run in the conversion workbook, this will again need to be changed manually (or scripted) according to the document 'Trout Brook RT Entrnace 4 2018 Designations'
entrance_class_join <- entrance_classifications %>%
  mutate(Entrance = as.character(Entrance)) %>%
  mutate(Start_time = format(strptime(substr(as.POSIXct(sprintf("%04.f", as.numeric(Start_time)), format="%H%M"), 12, 16), '%H:%M'), '%I:%M %p'))

step_III <- left_join(step_IIa, entrance_class_join, by = c('Park_id' = 'Park_trail', 'Entrance', 'Start_time', 'Day_type'))

# Check that there's no missing data for class entrances, as all entrances should receive a classification.  If there are observations missing a classification, this is most likely due to incorrect data entry (conversion from sample class workbook to index has been audited multiple ways and de-bugged; therefore it is shown to have captured 100% of the data in the sample class workbooks).  Check if trail entrance exists; if trail entrance does not exist in Sample Class Workbook, check entrance descriptions.  If it doesn't exist in entrance descriptions, most likely this is data error.  Contact IA coordinator for corrections.

missing_class_assignment <- step_III %>%
  filter(is.na(Usage_class)) %>%
  group_by(Park_trail, Agency_name, Start_time, Entrance) %>%
  count()

step_III_dups <- step_III %>%
  group_by(Park_id, Date, Start_time, Entrance, ResponseId) %>%
  count() %>%
  filter(n > 1)

step_III_dups_full <- inner_join(step_III, step_III_dups, by = c('Park_id', 'Date', 'Start_time', 'Entrance', 'ResponseId'))

Step 4: Add persons per vehicle and persons per bus multipliers, defined by the 2016 Visitor Study and 1998-1999 Visitor Study, respectively

step_IV <- left_join(step_III, bus_vehicle_multipliers, by = c('Park_id' = 'Park'))

Code below checks which parks did not have multiplier

step_IV_unjoined <- step_IV %>%
   filter(is.na(ppb_multiplier) | is.na(ppv_multiplier)) %>%
   select(Park_id) %>%
   unique()

## These parks are new and therefore don't have assigned multipliers; however they have zero vehicles or buses entering the park

# Check full number of observations missing data
step_IV_unjoined_full <- step_IV %>%
  filter(is.na(ppb_multiplier) | is.na(ppv_multiplier))

# Check number of buses and vehicles that came into the new parks
step_IV_unjoined_bv <- step_IV %>%
  filter(Park_id == 133 | Park_id == 140) %>%
  select(Vehicles, Buses) %>%
  filter(Vehicles != 0 | Buses != 0)


## Change missing ppv and ppb multipliers to 1
step_IV_tidy <- step_IV %>%
  mutate(ppb_multiplier = ifelse(is.na(ppb_multiplier), 1, ppb_multiplier),
         ppv_multiplier = ifelse(is.na(ppv_multiplier), 1, ppv_multiplier))

Step 5: Calculate expansion factors & attribute to parks

expansion_factors <- entrance_classifications %>%
  group_by(Park_trail, Usage_class, Day_type) %>%
  count() %>%
  rename(Expansion_factor = n)

# Join multipliers to parks to use to calculate total number of visitors in sample

step_V <- left_join(step_IV_tidy, expansion_factors, by = c('Park_id' = 'Park_trail', 'Usage_class' = 'Usage_class', 'Day_type' = 'Day_type'))


# Check that observations are not missing an expansion factor assignment.  No observation should be missing an assignment (i.e. this code should return a tibble with 0 observations).  If observations are missing an assignment, follow the same steps outlined in the missing assignment check of Step 3 of this script.
step_V_unjoined <- step_V %>%
  filter(is.na(Expansion_factor))

# For purposes of preliminary check of numbers, use only those counts that have an expansion factor assigned (un-assigned factors may be because parks are new or because of incorrect data entry)
step_V_not_na <- step_V %>%
  filter(!is.na(Expansion_factor))

Step 6: Compute expanded total of visitors (total for entire season, based on sample) at entrance and by date and time of day

# Compute expanded total of visitors (total for entire season, based on sample)

step_VI <- step_V_not_na %>%
  mutate(Expanded_total = (ppv_multiplier*as.numeric(Vehicles) + as.numeric(Bikers) + as.numeric(Peds_skaters) + ppb_multiplier*as.numeric(Buses) + as.numeric(Horse_riders) + as.numeric(Boaters)) * as.numeric(Total_days_by_type) * as.numeric(Expansion_factor))

Step 7: Add to previous three years of data

setwd('Datasets/2014-2017')
x2014_2017 <- read_csv('2014-2017 Data.csv')

x2015_2017 <- x2014_2017 %>%
  mutate(Least_recent_year = min(Year)) %>%
  filter(Year != Least_recent_year) %>% # Filter down to three most recent years of data
  select(-Least_recent_year) %>%
  separate(gooddate, into = c('Month', 'Day'), sep = 1) %>%
  mutate(Year_for_date = Year) %>%
  unite(Date, Month, Day, Year_for_date, sep = '-') %>%
  mutate(Date = as_date(mdy(Date))) %>%
  rename(Agency_id = agency,
         Park_id = Park,
         Bikers = Biker,
         Peds_skaters = Pedskat,
         Boaters = Boat,
         Vehicles = Vehicle,
         Buses = Bus,
         Horse_riders = Horse,
         Usage_class = samclass,
         Total_visitors = totalall,
         Day_type = daytype,
         Start_time = orgstart,
         Expansion_factor = expfact,
         ppv_multiplier = new2008ppv,
         ppb_multiplier = newppb,
         Total_days_by_type = `day#s`,
         Half_of_summer = half,
         Expanded_total = totalsum) %>%
  select(-carpeop, -buspeop) %>%
  mutate(Start_time = str_replace(Start_time, "10000", "1000"),
         Start_time = as.numeric(Start_time)) %>%
  mutate(Usage_class = ifelse(Usage_class == 1, 'High',
                              ifelse(Usage_class == 2, 'Medium', 'Low'))) %>%
  mutate(Start_time = format(strptime(substr(as.POSIXct(sprintf("%04.f", as.numeric(Start_time)), format="%H%M"), 12, 16), '%H:%M'), '%I:%M %p')) %>%
  mutate(Day_type = ifelse(Day_type == 1, 'Weekend_holiday', 'Weekday'))


x2018 <- step_VI %>%
  mutate(Year = 2018) %>%
  mutate(Entrance = as.numeric(Entrance),
         Bikers = as.numeric(Bikers),
         Peds_skaters = as.numeric(Peds_skaters),
         Vehicles = as.numeric(Vehicles),
         Buses = as.numeric(Buses),
         Boaters = as.numeric(Boaters),
         Horse_riders = as.numeric(Horse_riders))
  

x2015_2018 <- bind_rows(x2015_2017, x2018)

Step 8: Compute averages by Park, Trail Entrance/Timeslot Classification, and Day type (Weekend/holiday or Weekday)

x2015_2018_avg <- x2015_2018 %>%
  mutate(Park_id = ifelse(Park_id == 121, 116, Park_id)) %>% # Not sure why, but in past years SW RT was coded as 121.  In all documentation, however, it's coded as 116
  unique() %>%
  group_by(Park_id, Usage_class, Day_type) %>%
  mutate(Average_expanded_total = mean(Expanded_total)) %>%
  ungroup() %>%
  select(Park_id, Usage_class, Day_type, Average_expanded_total) %>%
  unique() %>%
  select(-Usage_class, -Day_type) %>%
  group_by(Park_id) %>%
  mutate(Estimate_18 = sum(Average_expanded_total)) %>%
  select(-Average_expanded_total) %>%
  unique() %>%
  ungroup() %>%
  mutate(Park_id = as.numeric(Park_id)) %>%
  arrange(Park_id)

Step 8a: Compare to previous year’s estimates; note that large differences do not inherently mean incorrect estimates

setwd('Datasets/2014-2017')
estimates_2017 <- read_csv('2017 Use Estimates.csv')

x2017_tidy <- estimates_2017 %>%
  rename(Park_trail = `Park or Trail`,
         Estimate_17 = `Total by class and day type`) %>%
  separate(Park_trail, into = c('Park_id', 'Trail_name'), sep = '=') %>%
  mutate(Park_id = trimws(Park_id),
         Trail_name = trimws(Trail_name)) %>%
  mutate(Park_id = ifelse(Park_id == 121, 116, Park_id)) %>%
  mutate(Park_id = as.numeric(Park_id))

estimates_17_18 <- left_join(x2015_2018_avg, x2017_tidy, by = c('Park_id'))

estimate_comparison <- estimates_17_18 %>%
  mutate(Estimate_17 = ifelse(is.na(Estimate_17), 0, Estimate_17)) %>%
  mutate(Estimate_gross_difference = Estimate_18-Estimate_17) %>%
  mutate(Percent_of_18 = (Estimate_18/sum(Estimate_18))*100,
         Percent_of_17 = (Estimate_17/sum(Estimate_17))*100,
         Estimate_percent_of_whole_difference = Percent_of_18-Percent_of_17) %>%
  mutate(Percent_difference_from_17_18_by_trail = 100-((Estimate_18/Estimate_17)*100))

# Percent difference in Trail or Park's Share of All Regional Visits between 2017 and 2018
ggplot(estimate_comparison, aes(Trail_name, Estimate_percent_of_whole_difference)) +
  geom_bar(stat = "identity") +
  coord_flip() +
  labs(title = "Percent Difference in Trail or Park's Share of all Regional Visits from 2017 to 2018")

# Percent Difference from 2017 Estimate to 2018 Estimate
ggplot(estimate_comparison, aes(Trail_name, Percent_difference_from_17_18_by_trail)) +
  geom_bar(stat = "identity") +
  coord_flip() +
  labs(title = "Percent Difference in Trail or Park's 2017 Estimate to Its 2018 Estimate")

# Gross Difference from 2017 Estimate to 2018 Estimate
ggplot(estimate_comparison, aes(Trail_name, Estimate_gross_difference)) +
  geom_bar(stat = "identity") +
  coord_flip() +
  labs(title = "Gross Difference from 2017 Estimate to 2018 Estimate")

estimate_comparison %>%
  gather(Estimate_17, Estimate_18, key = 'Estimate_year', value = 'Estimate') %>%
  ggplot(aes(Trail_name, Estimate, fill = Estimate_year)) +
  geom_bar(stat = 'identity', position = 'dodge') + coord_flip()

---

Following is an exploration of outliers

# Plot 2018 counts as striplot; use plotly to mouse over specific observations that appear to be outliers
stripchart_18 <- x2018 %>%
  mutate(Total_visitors = (ppv_multiplier*as.numeric(Vehicles) + as.numeric(Bikers) + as.numeric(Peds_skaters) + ppb_multiplier*as.numeric(Buses) + as.numeric(Horse_riders) + as.numeric(Boaters))) %>%
  ggplot(aes(Park_id, Total_visitors, color = Agency_id)) +
  geom_point(alpha = 0.6, size = 2) +
  coord_flip() +
  scale_color_distiller(palette = "RdPu") +
  labs(title = '2018 Raw Counts',
       y = 'Visitors Counted') +
  theme(axis.text.y = element_text(size = 8))

ggplotly(stripchart_18)

# Plot 2018 counts as boxplot; use plotly to zoom & mouse over specific observations that appear to be outliers
boxplot_18 <- x2018 %>%
  mutate(Total_visitors = (ppv_multiplier*as.numeric(Vehicles) + as.numeric(Bikers) + as.numeric(Peds_skaters) + ppb_multiplier*as.numeric(Buses) + as.numeric(Horse_riders) + as.numeric(Boaters))) %>%
  ggplot(aes(Park_trail, Total_visitors)) +
  geom_boxplot(alpha = 0.6, size = .5) +
  coord_flip() +
  labs(title = '2018 ',
       y = 'Counted Visitors') +
  theme(axis.text.y = element_text(size = 8))

ggplotly(boxplot_18)

# Plot 2018 distribution curve of Silverwood SRF by day type
 x2018 %>%
  filter(Park_trail == 'Silverwood Special Recreation Feature') %>%
  mutate(Total_visitors = (ppv_multiplier*as.numeric(Vehicles) + as.numeric(Bikers) + as.numeric(Peds_skaters) + ppb_multiplier*as.numeric(Buses) + as.numeric(Horse_riders) + as.numeric(Boaters))) %>%
  ggplot(aes(Total_visitors, color = Day_type, fill = Day_type)) +
  geom_density(alpha = 0.1) +
  labs(title = '2018 ',
       x = 'Counted visitors',
       y = 'Observations') +
  theme(axis.text.y = element_text(size = 8))

x2018 %>%
  filter(Park_trail == 'Silverwood Special Recreation Feature') %>%
  mutate(Total_visitors = (ppv_multiplier*as.numeric(Vehicles) + as.numeric(Bikers) + as.numeric(Peds_skaters) + ppb_multiplier*as.numeric(Buses) + as.numeric(Horse_riders) + as.numeric(Boaters))) %>%
  ggplot(aes(Total_visitors, fill = Day_type)) +
  geom_histogram(size = 1) +
  labs(title = '2018 ',
       x = 'Counted visitors',
       y = 'Observations') +
  theme(axis.text.y = element_text(size = 8))

# Plot 2018 counts for Silverwood SRF as density curve, no day type
 x2018 %>%
  filter(Park_trail == 'Silverwood Special Recreation Feature') %>%
  mutate(Total_visitors = (ppv_multiplier*as.numeric(Vehicles) + as.numeric(Bikers) + as.numeric(Peds_skaters) + ppb_multiplier*as.numeric(Buses) + as.numeric(Horse_riders) + as.numeric(Boaters))) %>%
  ggplot(aes(Total_visitors)) +
  geom_density() +
  labs(title = '2018 ',
       x = 'Counted visitors',
       y = 'Observations') +
  theme(axis.text.y = element_text(size = 8))

 # Plot all trails/parks density curves 
x2018 %>%
  mutate(Total_visitors = (ppv_multiplier*as.numeric(Vehicles) + as.numeric(Bikers) + as.numeric(Peds_skaters) + ppb_multiplier*as.numeric(Buses) + as.numeric(Horse_riders) + as.numeric(Boaters))) %>%
  ggplot(aes(Total_visitors)) +
  geom_density() +
  labs(title = '2018 ',
       x = 'Counted visitors',
       y = 'Observations') +
  theme(axis.text.y = element_text(size = 8)) +
  facet_wrap(~Park_trail)

ggsave('2018 Distribution Curves for Counted Visitors.jpg', width = 50, height = 30, limitsize = FALSE)